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Abstract 

This paper describes a package for calculations of expressions with 
Dirac matrixes. Advantages to existing similar packages are described. 
MatrixExp package is intended for simplification of complex expres- 
sions involving 7-matrixes, providing such tools as automatic Feyn- 
man parameterization, integration in d-dimensional space, sorting and 
grouping of results in a given order. Also, in comparison with existing 
similar package Tracer, presented package MatrixExp has more en- 
hanced input possibility. User-available functions of MatrixExp pack- 
age are described in detail. Also an example of calculation of Feynman 
diagram for process b — > s^g with application of functions of Matrix- 
Exp package is presented. 
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PROGRAM SUMMARY 

Keywords: Computer algebra, Particle Physics, Computing, Quantum Field The- 
ory, High energy physics. Calculation, Gamma-algebra, Dimensional Regularization 
Classification: 11.1 General, High Energy Physics and Computing, 4.4 Feynman 
diagrams 

Nature of problem: 

Feynman diagram calculation, simplification of expressions with 7-matrixes 

Restrictions: 

MatrixExp package works only with single line of expressions (G[ll,...]), in con- 
trary with Tracer package which works with multiple lines, i.e. the following is 
possible in Tracer, but not in MatrixExp : G[11,...]**G[12,...]**G[13,...]..., which 
will return the result of G[11,...]**G[11,...]**G[11,. ..]... 

Running time: 

Seconds for expressions with several different 7-matrixes on Pentium IV 1.8GHz 
and of the order of a minute on Pentium H 233MHz. Calculation times rise with 
the number of matrixes. 
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LONG WRITE-UP 



1 Introduction 

Modern calculations of radiative corrections in Relativistic Quantum Field 
Theory, in particular within the limits of Standard Model, lead to calculations 
of one-, two- and higher-loop Feynman diagrams. Considering, that calcula- 
tions of two-loop and furthermore multi-loop diagrams are lab or- consuming 
enough and volumetric, and also that the number of possible diagrams to be 
evaluated is great, now are widely applied various packages of symbolic ma- 
nipulations. In particular. Tracer package [T] working in MATHEMATICA 
[2j environment is known, allowing to operate with expressions of 7-matrixes 
in (i- dimensional space. Following the developers of Tracer package we chose 
MATHEMATICA environment for the MatrixExp package, since MATHE- 
MATICA is widely used by researchers and is easily programmable, and also, 
which is more important, unlike similar package MAPLE MATHEMAT- 
ICA allows to program in "rule-based" style: a special style of programming 
convenient for implementation (application) of mathematical knowledge into 
the program (package). From the technical point of view, during calculations 
of radiative corrections it is necessary to get rid from arising intermediate 
singularities, inherent to quantum corrections, usually by the method of di- 
mensional regularization; hence, it is necessary to perform the calculations 
in d-dimensional space. The MatrixExp package uses the so-called Naive Di- 
mensional Regularization scheme (NDR). As it is known, the given scheme 
leads to some algebraic inconsistencies in (i- dimensional space jU |3] . How- 
ever, as it has been shown in Ref. jH], in many calculations which don't 
involve traces with 75 matrix, NDR scheme gives correct result. A problem 
connected with definition of 75 matrix in dimensional space is solved in the 
t'Hooft-Veltman regularization scheme [71 |H]. In future, this scheme will also 
be included in the MatrixExp package. 

2 Tracer package 

Tracer package is written basically in "rule-based" style which is presented 
in pp. We shall briefly describe this style since our developed program also 
is written in this style. Unlike the procedural style of programming, where 
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the program is executed consistently from the first command up to the fi- 
nal, realizing the algorithm for solving the task in view, "rule-based" style 
is setting up a collection of transformation rules which just represent math- 
ematical rules of transformations and properties of the given objects. The 
decision of a task in view in "rule-based" style is reduced to the following: 

1. Setting up of transformation rules (which are applied when expression 
matches the pattern on the left-hand side (Ihs) of the rule, and return 
the right-hand side (rhs) of the rule) 

2. Use of basic rules buih-in in MATHEMATICA 

3. Consecutive application of corresponding rules (automatically carried 
out by MATHEMATICA) 

4. Last received expression which does not match to any pattern is the 
received result. 

Schematically the solution of a problem in "rule-based" style is presented 
on Fig. 1. 

The solution of a problem in this style greatly simplifies the programming 
and does clearer the program code since allows to not supervise all stream of 
the information during the solution of a problem (at calculation of diagrams 
very bulky and complex expressions arise), but by means of setting up of all 
transformation rules (which are the initial mathematical rules, i.e. are clear 
and natural to the user) directly solves the problem in view (in our case - will 
transform initial expression to a convenient form, carrying out all necessary 
transformations). To illustrate the above mentioned, codes for differentiation 
in procedural style and in "rule-based" style are presented in Table ^ 

Note, that if we try to define differentiation for all^ other functions as 
well, then the code in procedural style would get very complicated, including 
special recursive codes for mixed functions, code for parsing expression of 
functions etc., and in the case of "rule-based" programming, one has just 
to add all the definitions, rules for differentiation of a sum and product of 
functions and rules for differentiation of f{g{x)). 

Tracer package realizes 7-algebra in "rule-based" style and calculates a 
trace of strings of 7-matrixes [T] . Tracer accepts following types of data in a 
corresponding format: 

^e.g. for 5 or 10 more functons 
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Transformation rules, specified by 
the programmer 



Basic transformation rules (implemented in 
MATHEMATICA) 



all transformation rules 



I 



expression to calculate (input) 




Figure 1: Scheme of problem solution in a "rule-based" style. 
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a. Procedural style 

diff[y_,x_] : = 
Block[{n}, 

If[TrueQ[Length[y]==2 && y[[0]]==Power && y[[l]]==x], 

n=y[[2]]; n*x' (n-1), 
If[TrueQ[y==x], 1, 

If[TrueQ[Length[y==l] && y[[0]]==Log && y[[l]]==x], 1/x]]]; 

b. "rule- based" style 

diff[x_' n_, x_] := n*x" (n-1); 
diff[Log[x_],x_]:= 1/x; 



Table 1: Codes for derivative calculation (cases of power function and loga- 
rithmic function) in procedural (a) and "rule-based" (b) styles. 

1. a momentum or a sum of momenta p, p — q, {p + kl — k2), which 
correspond to p, p — q and p + ki — k2 

2. an open index of a gamma matrix, which is denoted by a symbol in 
curly brackets {alpha} corresponding to 7o, 

3. a linear combination of the unit matrix in the gamma algebra U and 
the matrix 75 denoted by G5, with scalar coefficients 

4. a linear combination of a momentum p and a mass term mU. 

For convenience, we use the same input-format for MatrixExp package^, 
for users already familiarized using Tracer package. Also, this way one may 
directly use results achieved by Tracer as input for MatrixExp and vice versa. 

''except points (3) and (4). There are special symbols for (1 ±75)72 defined as R and L 
correspondingly. 75 also may be used in the same way, i.e. G5. For point (4) see sections 
3 and 4. 
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3 MatrixExp Package 



When calculating diagrams, usually so-called Feynman parameterization jH] 
is applied, therefore, in initial expression instead of momenta p arises their 
linear combinations, for example p + ui* k + U2* r... and so the direct appli- 
cation of Tracer package is already impossible, i.e. it is necessary to segment 
(partition) the initial expression to parts, not containing such combinations 
of momenta {p,k,r...) and scalars (mi,M2...). 

Developed new package MatrixExp presented in this paper, allows enter- 
ing data without additional transformations. 

Input data for calculation by Tracer and MatrixExp packages correspond- 
ingly (b) and (c), is presented on Table [21 

Thus, MatrixExp package allows working with complex expressions, and 
there is no need for additional manual transformations of initial input expres- 
sion. Besides, since it is very often applied Feynman parameterization, with 
the subsequent d- dimensional integration, this opportunity is also included 
in the MatrixExp package, i.e. it is possible to set both automatic Feynman 
parameterization, and automatic integration. The package is intended not 
only for calculation of traces of matrixes, but also for calculations and sim- 
plifications (sorting and grouping) of complex expressions received as a result 
of previous calculations or as output from other packages. During calcula- 
tions MatrixExp package operates with scalar products (additional conditions 
imposed on scalar products can be set via appropriate user-function of Ma- 
trixExp (see next section)), combines and recombines momenta according to 
p = '^aPa for performing calculations, integrations and grouping results (for 
example, Paghpiuir^ + 'Jahfi--- ^ q{ui{p.r) +p)k..., etc.). 

MatrixExp package includes such enhancements as definition of scalars, 
automatic sorting by the given mask (order), automatic Feynman param- 
eterization (using previously declared scalar symbols), grouping of sorted 
results, automatic (i-dimensional integration over specified momenta, auto- 
matic rearrangement of "edge" -momenta {ph and Ps) and their replacement 
by corresponding masses according to Dirac equation {{s\ps = ms{s\ and 
Pb\b) = mb|6)), operations with scalar products, etc. 
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(a) Expression to be calculated 

{s\Ps: la:r,p-mb,k + f,-fa\b), T - UiP + U2Ps 

(b) Input for Tracer package 

all = G[ll, ps, {alpha}, r, p-U mb, k+r, {alpha}]; 
al2 = -ul*G[ll, ps, {alpha}, k, p-U mb, k+r, {alpha}]; 
al3 = u2*G[ll, ps, {alpha}, ps, p-U mb, k+r, {alpha}]; 
a21 = -ul*G[ll, ps, {alpha}, r, p-U mb, k-p, {alpha}]; 
a22 = ul^2*G[ll, ps, {alpha}, k, p-U mb, k-p, {alpha}]; 
a23 = -ul*u2*G[ll, ps, {alpha}, ps, p-U mb, k-p, {alpha}]; 
a31 = u2*G[ll, ps, {alpha}, r, p-U mb, k+ps, {alpha}]; 
a32 = -u2*ul*G[ll, ps, {alpha}, k, p-U mb, k+ps, {alpha}]; 
a33 = u2" 2*G[11, ps, {alpha}, ps, p-U mb, k+ps, {alpha}]; 
result = all+al2+al3+a21+a22+a23+a31+a32+a33; 

(c) Input for MatrixExp package 

Vsetscalars[ul, u2, mb]; 

a = G[ll, ps, {alpha}, r-ul*p+u2*ps, p-mb, k+r-ul*p+u2*ps, {alpha}]; 
result = Vcalc[a]; 



Table 2: Galculation of expression (a) using Tracer package (b) (as we see, 
one needs to segment the initial expression into parts, and if the expression 
is more complex, then the calculation becomes very complicated) and using 
MatrixExp package (c) (Function Vsetscalars[] declares -ul, u2 and mh as 
scalars and then function Vcalc automatically performs all corresponding 
operations). 

4 User functions of MatrixExp package 

The detailed description of MatrixExp functions and their syntax is included 
in the package and is accessible via standard method in MATHEMATICA, 
i.e. ' ?function_name '. Call ' Vhelp[] ' to receive the list of all functions 
defined in the package. 

Here we present basic functions and their purpose. 

Vsetscalars[ul, u2...] - declares symbols ul,u2... as scalars. By default 
the package declares as scalars the following symbols: ul...u9 since they 
are necessary to perform Feynman parameterization. It is possible to rede- 
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fine scalars calling function Vsetscalars [...], or to add other scalars calling 
function Vaddscalars [...]. 

Vsetintvar[r] - declares the symbol r as integration variable over which 
Feynman parameterization and automatic ci-dimensional integration will be 
performed. 

Vsetsortlist[pl, {mul}, {mu2}, p2, p3...] - sorts momenta (pi, p2, p3...) 
and 7-matrixes ({mul}, {mu2}...) in the designated sequence. 

Vcalc[exp, options] - carries out calculations on expression exp and 
groups (sorts) result. The following may be specified as control parameters 
(options): integration - On/Off (by default it is switched off), a symbol of 
integration (integrating — True/False, integrating — > r) (the latter also 
implies integration True), Feynman parameterization On/Off (by de- 
fault it is switched off) (dofeynman — > True/False), sorting On/Off, sort- 
ing order (sorting — > True/False, sorting —>■ {p, {mu}, k, {nu}, {sigma}}) 
(the latter also implies sorting True), rearrangement of "edge" -momenta 
On/Off (by default it is switched off) (pps True/False), calculation of a 
trace On/Off. (by default it is switched off) (spur —>■ True/False)...., etc. 

Vfeynman[x] - performs Feynman parameterization of expression x. 
Prior to calling of this function it is necessary to define the variable of integra- 
tion by calling Vsetintvar[r], with respect to which function Vfeynman[x] 
should make parameterization. Gives the output in the following form: {K, 
scalar'i scalar 2 ...scalar^^ , int_varJ}, power, delta, N}, where - is the fac- 
tor arising in a result of parameterization, scalar^ , scalar'^2 i scalar^^ - are 
additional multipliers, if any arise, intjuarS) - is the displacement of the in- 
tegration variable, i.e. in initial expression the integration variable (let it be 
X for example) should be replaced by X — int_varJ}, power - is the power 
of the resulting denominator after parameterization, delta - is the parame- 
ter of integration, iV - is not factored part, i.e. multipliers and factors not 
dependent on the integration variable, which have not been parameter izated 

m 

1 r duldu2...dun uV^...un'^'"'6{l — ul — ... — un) 

'dTD2D3D¥ ~ J [(int_var - int_var_0)2 + delta]P°*'^'- ^ ^ 

When specifying Feynman parameterization option in function Vcalc, all 
steps and substitutions are carried out automatically (see expressions (|Z1E1))- 

Vsetrules[{pl, p2, cl}, {pi, c2}...] - defines rules of scalar product, cor- 
responding to - ((pl,p2) = cl, {pi, pi) = pl^ = c2). 
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5 Calculation of an expression by means of 
MatrixExp package 

Below we present an example of calculation with use of MatrixExp package. 
It is calculated a Feynman diagram for process b —>■ s^g ^U] in case of local 
operators Oi,02, having the following form: Oi = {sL'~^^T°-CL){cL'-)^T°'bL), 
O2 = {sL'~i^j.CL){cL'^^hL). On Fig. 2 we present the two Feynman diagrams 
sum of which in |I0; is denoted as Jafs- Here q designates photon momentum, 
r - gluon momentum. 




Pa a p 

a) b) 



Figure 2: Bremsstrahlung diagrams induced by Oi and O2. 



- lfM,k + q + mc,ja,k + mc,-ff},k-f + mc,jf,, 

Ml ~ sl— — -ol, [2) 

[(k + qY — T^c] [(k — rY — Tnf\ [k? — m^] 



M2 ~ Sl— — ; — -Ol, (3) 

[{k + f)^ — m^] [(k — qY ~ m?^ [k"^ — m^] 

J^p = M1 + M2. (4) 

On start-up MatrixExp package by default defines some parameters, in par- 
ticular variables ul...u9 are defined as scalars, automatic sorting is switched 
On, automatic integration is switched Off. 
First of all we add mc as scalar 

Vaddscalarsfmc] (5) 
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Next we define the rules for scalar products = 0, = 

Vsetrules[{q,0},{r,0}] (6) 
Finally we call function Vcalc for calculation of Mi and M2 accordingly 

G[ll, R, {mu}, k + q + mc, {al}, k + mc, {bet}, k — r + mc, {mu}, L] 



Ml = Vcalc 



((k — r)^ — mc^)((k + q)^ — mc^)(k^ — mc^) ' 
integrating k, dofeynman True, sorting {{al}, {bet}, r, q}] (7) 

M2 = Vcalcf ^^^"*^' ^' k + r + mc, {bet}, k + mc, {al}, k - q + mc, {mu}, L] 

((k + r)2 — mc2)((k — q)^ — mc'^){}i? — mc^) 
integrating k, dofeynman True, sorting {{al}, {bet}, r, q}] (8) 

where the option dofeynman True turns On Feynman parameteriza- 
tion, and the option integrating k defines the integration variable (and 
turns On integration). Received Ml and M2 are grouped in sorted re- 
sulting matrix expressions (e.g. G[ll, {al}, {bet}, q], G[ll, {al}, {bet}, r], 
G[ll, {al}, {bet}], G[ll, {bet}, r, q]...) and include scalars m1, u2 and m3, 
which arose due to Feynman parameterization, (actually one of the scalars 
will be missing (in this case - ul) since during Feynman parameterization, 
MatrixExp package automatically carries out integration (/ X ■ 5(1 — ul — 
u2 — u3)dul), replacing ul by 1 — m2 — u3). After that we consider the sum 
Jal3 = M1 + M2 

Jab = Vcalc[Ml + M2, sorting {{al}, {bet}, r, q}] (9) 

or, using the built in function of MATHEMATICA (since Ml and M2 are 
already sorted in the identical order) 

Jab = Simplify [Ml + M2] (10) 

In the further we need to simplify expression (jHl ITn|l , using transversality 
of photon and gluon. Ward identities, etc., and after representing through 
E[a,(3,r] we receive the form presented in |lUj . 
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TEST RUN OUTPUT 



In the package d = 4: — 2e denotes the space-dimension, so avoid its usage 
as momenta.^ 

In[361] := Vcalc[G[ll, a, {mu}, b, {mu}, a]] 
Out[361] = (-2 + d)G[ll,b]a.a-2(-2 + d)G[ll,a]a.b (11) 
Timing := 0.05 seconds 

To be able to use scalars in expressions we must declare scalars, otherwise 
by default they will be treated as momenta 

In[362] := Vaddscalairs[ma, mb, mc] 

Out [362] = {ul, u2, u3, u4, u5, u6, u7, u8, u9, ma, mb, mc} (12) 
Timing := 0.00 seconds 

We can use either global sorting (via Vsorton[True] and Vsetsortlist[a, b, c], 
which will be automatically the default sorting for all expressions) or local 
sorting 

In[363] := rO = Vcalc[G[ll, R, a + ma, b, c + mc, L]] 
Out [363] = mamcG[ll,b,L] + G[ll,a,b,c,L] (13) 
Timing := 0.00 seconds 



In[364] := rl = Vcalc[G[ll, R, a + ma, b, c + mc, L], 

sorting — > {c, b, a}] 
Out[364] = -G[ll,c,b,a,L] + 2G[ll,c,L]a.b + 

G[ll, b, L] (mamc - 2 a.c) + 2 G[ll, a, L] b.c (14) 
Timing := 0.02 seconds 



In[365] := r2 = Vcalc[rl, sorting — > {b, c, a}] 

Out [365] = G[ll, b, c, a, L] + G[ll, c, L] a.b+ 

G[ll,b,L] (mamc-2a.c) (15) 

Timing := 0.00 seconds 

^Thc last row "Timing: " shows the tested time of the calculation. Results for Pentium 
IV 2.8GHz. 
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In[366] := r3 = Vcalc[r2, sorting {a, b, c}] 

Out [366] = mamcG[ll,b,L] + G[ll,a,b,c,L] (16) 
Timing := 0.02 seconds 

Here L and R denote (1 — 75)/2 and (1 + 75)/2 respectively. As we see, the 
Output (including numerical coefficients and scalars) may be used as Input 
without any additional transformations (e.g. rl, which is Out [364] is used as 
input in (fT^ ). 

To use integration and Feynman parameterization, lets define some rules 
for scalar product (we use arbitrary expressions and definitions here, just as 
mathematical expressions) 

In[367] := Vsetrules[{a, ma^}, {b, mb^}, {a, b, ma * mb/2}] 

Out[367] = {{a, ma^}, {b, mb^}, {a, b, ^^1^}} (17) 
Timing := 0.02 seconds 



ln[368] .- res - Vcalc[ ((j, _ _ ^a2)((k - b)^ - mb^) ' 
integrating — k, dofeynman True] 

Out[368]= -((-2 + d)mb(47r)-'+^P^ 

/ \ -l+eps 

\ -ma2(-l + u2)2 + mamb (-1 + u2) - mb2u22y 
(-(-2 + eps) (-ma^(-l + u2)^ + mamb (-1 + u2) u2- 
mb^ u22)) + (-1 + eps) (ma^ (-1 + u2)^+ 
mamb (-1 + u2) u2) - mb^u22) + (-1 + eps) 
(ma^(-l + u2)2 - mamb (-1 + u2) u2 + mb^ u2^)) 
G[ll, a, L] Gamma[-1 + eps])/(-ma^ (-1 + u2f+ 
mamb(-l + u2)u2-mb2u22) (18) 
Timing := 0.13 seconds 

We see, that all the substitutions after parameterization are done, a Dirac- 
delta integration over parameter ul is done, integration over k is done, and 
we have a final result. 
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Here we calculate the spur in two ways. First we calculate the expression 
rO, followed by calculation of the spur of the result, and then we directly 
calculate the spur of the initial expression rO. 



In[369] := rO = G[ll, a + ma, b + mb, a + ma, c + mc, f + ma, L]; 
rr = Vcalc[rO] 
Vcalc[rr, spur True] 
Vcalc[rO, spur True] 

Out[370] = 3ma^mbmcG[ll,L] + 3ma2mbmcG[ll,a,L] + 
3 ma'^ mb G [11 , c, L] +3 ma^ mb mc G [11 , f , L] + 
3 ma^ mb G[ll, a, c, L] + 3 mamb mc G[ll, a, f, L] + 
3 ma^ mb G[ll, c, f , L] + 3 mamb G[ll, a, c, f, L] 

Out[371] = 6mamb (maa.c + mca.f + ma (mamc + c.f)) 

Out [372] = 6 mamb (maa.c + mca.f + ma (mamc + c.f)) (19) 

Timing := 0.09 seconds 

And finally here are the timings for calculations of ((THSl): 1-86, 1.64 and 
6.00 seconds respectively. 
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